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Abstract 

The choice of a point set, to be used in numerical integration, determines, to a 
large extent, the error estimate of the integral. Point sets can be characterized 
by their discrepancy, which is a measure of its non-uniformity. Point sets with 
a discrepancy that is low with respect to the expected value for truly random 
point sets, are generally thought to be desirable. A low value of the discrepancy 
implies a negative correlation between the points, which may be usefully employed 
to improve the error estimate of a numerical integral based on the point set. 
We apply the formalism developed in a previous publication to compute this 
correlation for one-dimensional point sets, using a few different definitions of 
discrepancy. 
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1 Introduction 



In a previous publication 0, we have discussed the problem of computing a nu- 
merical integral estimate and its error, using a given set consisting of N points Xk, 
k = 1, 2, . . . , N in the multi-dimensional hypercube K. As is well known, in the 
Monte Carlo approach, where the points are chosen iid random with the uniform 
distribution, the expected error decreases as l/y/N. This can be improved upon 
by applying quasi-random point sets, that are, owing to some more or less so- 
phisticated generation algorithm, more uniformly distributed than random ones. 
This is embodied in the notion of discrepancy D^, a measure of the deviation 
from uniformity. A very uniformly distributed point set will have a relatively 
low value s of this discrepancy, whereas very non-uniform ones will have a high 
value of s. Obviously, different definitions of discrepancy are possible in addi- 
tion to the best-known so-called star- discrepancy D* N . Good uniformity implies, 
however, that the points are not independent, and in particular the distributions 
P\(x\) of a single point x\, and P 2 (xi,x 2 ), of a pair of points Xi,x 2 , will deviate 
from unity0: we can define 

PiM = l-^F 1 (s;x 1 ) , 

P 2 (x 1 ,x 2 ) = 1 - — F 2 (s;x 1 ,x 2 ) . (1) 

As discussed in [[jj, in order to talk sensibly about such distributions, we must 
have an underlying definition of an ensemble of point sets. In the truly random, 
Monte Carlo, case this is the ensemble of all sets of N points with the Cartesian 
product of N uniform one-point distributions: in the quasi-random case, we 
restrict this ensemble to all point sets having a fixed value s of the discrepancy. 
For truly random points, s is then a random variable, with a probability density 
H (s). In we have developed a formalism that allows us to define a useful 
form for Djy, and from that to compute a 1/N expansion for H (s), Fi(s; Xi) and 
F 2 (s; Xi, x 2 ). Let f(x) be a quadratically integrable function on K. For a general 
integration error 
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we can then show that, averaged over the restricted ensemble (under our definition 
of discrepancy), the integral is unbiased, i.e. (77) = 0, and 



N( V 2 ) = j dx(f(x)) 2 - Udxf(x) 

K XK j 

1 - -— ) J dxidx 2 f(x 1 )f(x 2 )F 2 (s; x 1: x 2 ) , (3) 



3 These are the distributions obtained by integrating over the remaining points. 
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so that we may expect an improved error if F 2 is positive whenever x\ and x 2 
are close together. For truly random points, the term with F 2 drops out, and we 
recover the standard Monte Carlo result. 



The treatment in JIJ has been purely formal; if it is to be any good, we are 
obliged to show how it works out in practice. The present paper is the first step 
in that direction: we shall apply the formalism to a few particular definitions of 
Djy in the case of one-dimensional point sets and integrals. Of course this is more 
or less academic since one-dimensional integrals can, in general, be approximated 
to extreme accuracy by other, non-stochastic methods. But, by keeping the 
extension to more dimensions in mind, we may hope to develop useful results and 
insights. The lay-out of this paper is as follows. In section 2, we briefly rehash 
the relevant results from 0. In section 3, we present results for a very simple, 
and not altogether too useful, discrepancy. In section 4, we proceed to a more 
realistic definition of discrepancy, which will lead to a surprising connection with 
well-known results for the star-discrepancy. We conclude with some preliminary 
remarks on the extension to the case of more- dimensional point sets and integrals. 



2 The formalism 

Given a point set of iV points xi, x 2 , ■ ■ ■ , x^ in the unit interval [0, 1), we define 
its discrepancy as 

1 N 

D N = T7 a l U n (Xk)u n (xi) , (4) 
iV n>0 k,l=l 

where u 2n (x) = ^/2cos(2irnx), u 2n -i(x) = \/2sin(27rna;), and the numbers a n 
depend on the class of integrands that our particular integrand is thought to be 
a member of. In the present model of integration, a typical integrand f{x) is of 
the form 

f(x) =V + J2 V nUn(x) , (5) 
n>0 

where the coefficients v n are distributed independently and normally around zero 
with standard deviation o n . We shall always assume translational invariance, 
that is, 02 n -i = cr 2ri . This implies that F 2 (s; xi, X2) only depends on x\ — x 2 , so 
that we may as well write it as F 2 (s; x% — x 2 ). Moreover, for the integrands to be 
quadratically integrable, the sum J2 n >o a n m ust converge. We define 



G (z) = exp(-W(l-2zo*)) 

\ Z n>0 / 



2za. 



) = E 1 _ 2 ; a2 M*K(0) • (6) 



n>0 
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Let now the value of Djv equal s for our point set. Then, we may compute H 
and F 2 to leading order in 1/N as follows: 



IOO 

H (s) = — [ dze- zs G (z) 



-IOO 
IDG 



R 2 (s;x) = -— I dz e zs (j)(z;x) G (z) , 

ZTTt J 

—ioo 

F2{S]X) = -^loW ' (?) 

where the integration contours run to the left of all singularities. Hence, both H Q 
and i?2 vanish trivially whenever s < 0. Owing to the translational invariance in 
our model, the function F\ vanishes to all orders in 1/N, making the integration 
unbiased. Moreover, the average value of Dn for truly random point sets is 

<Av> = E^ • (8) 

n>0 



Finally, from the definitions (|6|,0), we may immediately infer that 

which may serve as a check on the computations, and - since F 2 may be supposed 
to reach a maximum for X\ = x 2 , at least for small s - can yield an approximation 
to a lower bound on s, simply by putting F 2 (s; x, x) = N. 



3 A simple model 

The first, and simplest model, of integrands is obtained by taking 

< = ^ > n = l,2,...,M , (10) 

and zero for higher values of n. This means that our typical integrand has all 
Fourier components up to frequency M occurring with equal strength on the 
average, higher modes being absent. Note that this gives (s) = 1 for truly 
random point sets. In this case, we have 

z \ ~ M M M 

* n — • • 1V1 M-l-sM 



where the last result (for s > 0) follows when we close the integral by a contour 
to the right. Note that, for large M, Hq(s) approximates a Gaussian with unit 
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mean and width Xj^/M. The discrepancy is given by 

1 N 

D n = zl [Mx k - a;,) - 1] , 

k,l=l 

* < ^ . cos(M7rx)sin((M + 1)ttx) . . 

sm(7nr) 

A particularly interesting point set is the one with minimal discrepancy, namely 

x fc = {x + ^} , k = l,2,...,N , (13) 

where Xq is arbitrary, and {x} denotes x modi. The discrepancy depends, in this 
case, only on N and M: 

, M modiV M A , 

D N = l • (14) 

We see that, for M < N, the discrepancy actually vanishes. This is reasonable, 
since an integrand that only has Fourier modes of frequency iV — 1 or lower is 
indeed integrated with zero error by iV equidistant points. For large M, the 
value of D N lies generally quite close to one. The function F 2 (s; x) can easily be 
constructed in this case: since all the nonzero a n are equal, it must be proportional 
to £m{x) ~ I; an d from Eq.^ we can immediately find its normalization, to arrive 
at 

F 2 {s;x) = 2(1 - s) [£ M (z) - 1] • (15) 

In Fig. [l], we present £m( function of x for a few values of M. As M 

increases, approaches 1/2 + l/25({x}), where S(x) is the Dirac delta function. 
Finally, the requirement F 2 (s; 0) < N leads to s > 1 — N/ (2M), which is a quite 
useless bound if ^> M, but is reasonable if M ^> N. 

The main attraction of this model is its simplicity. The extension to more 
dimensions is straightforward. Also, as proven in the appendix to [|l|], one can 
indeed show that, relative to the expected error from classical Monte Carlo, the 
expected error obtained with a point set of discrepancy s is indeed reduced by a 
factor s. On the other hand, functions that have only a finite number of modes, 
and all equally strong at that, do not really form a good model for integrands 
likely to be encountered in practice. We therefore move, now, to a more realistic 
model. 



4 A non-simple model 

The next model we shall study is characterized by 

1 / X 

0"2n = - , 71=1,2,3,... . (16) 

n 
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Figure 1: The function £m(x) for various values of M. 

Therefore, all Fourier modes can occur in the integrand, the higher modes being 
typically smaller than the lowest ones. Note that, with respect to the summability 
requirement on the a n , this is not the slowest possible behaviour: we are, in 
principle, allowed to take a-in ~ n e_1 / 2 , with any positive e: but the above choice 
appears to be the most manageable. In fact, integrands with discontinuities 
typically have Fourier coefficients decreasing like 1/n, so that the above choice 
for a appears to be justified for integrands that may be discontinuous but in which 
quadratically integrable singularities have been mapped away. This is precisely 
the usual situation in, e.g., particle phenomenology. 

We may now reapply the formalism of section 2. In the first place, the dis- 
crepancy is given by 

1 N 

D N = T7 Z! Pl( x k-Xi) , 
iV k,l=l 

PiW = E ^cos(27rna;) = ^ [1 - 6{x}(l - {x})} , (17) 

n>0 71 3 

which follows trivially from the properties of Bernouilli polynomials. Some special 
cases are of interest. 

• The most non-uniform point set, in which all points coincide, has discrep- 
ancy D N = n 2 N/3. 

• For truly random points, the expected discrepancy is (D N ) = n 2 /3. 
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Generalized Richtmeyer sequences, defined by 

Xk = {x + ka} , (11 
with arbitrary xo and non- integer a, lead to a discrepancy 



TV 2 



D N = 

N 3N 



N-l 

N + 2 ( N ~ m )( x - 6{m«}(l - {ma})) 

m=l 



(19) 



which can conveniently be computed in time O (N). 

The most uniform point set, that of Eq.[T3, has a discrepancy 



The least and most uniform point sets are in fact, special cases of the generalized 
Richtmeyer sequence with a = and a = 1/N, respectively. A nice illustration 
is the following. We consider generalized Richtmeyer sequences, starting with 
a = 0. We then iterate a <— 1/(1 + a), which gives us a series of rational 
(Fibonacci) approximants to the golden ratio. For each a we then compute the 
discrepancy, for several values of N. The results are given in table 1. We see that 
the discrepancy decreases as the approximation improves: however, for every fixed 
iV there are some lower approximants that give a better discrepancy than some 
higher approximants: the phenomenon of irregularity of distribution. Similarly, 
with increasing N, the discrepancy decreases somewhat slower than with 1/N. 

The evaluation of H Q and F 2 is somewhat more complicated than in the pre- 
vious model. In the first place, we have 



^ V2 



G ^ = Vi-^r= . p^r ■ (21) 

n>o 71 ~ Zz sin V 271^2: 

The simplest way to compute H (s) is, then, to close the complex contour of the 
z integral to the right and contract it to run around all the poles at z — m 2 /2, 
m — 1, 2, . . .: this immediately gives us 

H (s) = £ (-) m ~ l m 2 e - sm2 l 2 . (22) 

m>0 

This result is, in some sense, surprising, since it is essentially identical to that of 
the star-discrepancy D* N : by comparison with the Kolmogorov-Smirnov law, we 
have 



Prob (D N <s) = Prob \D* N < yJs/(16N)J . (23) 

It is not completely clear why there should be such a simple relation between 
the two they are based on quite different models of the underlying class 
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iter. 


a 


JSI = 100 


W = 1000 


AT" 1 r\r\r\r\ 

N = 10000 


1 


1.00000000 


328.98681 


3289.86813 


32898.68134 


2 


.50000000 


82.24670 


822.46703 


8224.67033 


3 


/-I /-I /-I /-I /-I /-I /-I ^ 

.00000007 


36.58333 


365.54383 


3655.40933 


4 


.60000000 


13.15947 


131.59473 


1315.94725 


5 


.62500000 


5.18977 


51.40419 


514.04190 


6 


.61538462 


1.99806 


19.46995 


194.66713 


7 


.61904762 


.78420 


7.46330 


74.60073 


8 


.61764706 


.32284 


2.85192 


28.45961 


9 


.61818182 


.16748 


1.09343 


10.87618 


10 


.61797753 


.10234 


.41917 


4.15401 


11 


>~i -i n a r p p" z"" 1 

.61805556 


.09481 


.16257 


1.58743 


1 2 


61 802^7^ 


09296 


06628 


60668 


13 


.61803714 


.09298 


.02975 


.23233 


14 


.61803279 


.09287 


.01554 


.08931 


15 


.61803445 


.09290 


.00732 


.03473 


16 


.61803381 


.09288 


.00832 


.01404 


17 


.61803406 


.09289 


.00763 


.00604 


18 


.61803396 


.09289 


.00784 


.00303 


19 


.61803400 


.09289 


.00775 


.00206 


20 


.61803399 


.09289 


.00779 


.00146 


21 


.61803399 


.09289 


.00777 


.00155 


22 


.61803399 


.09289 


.00778 


.00149 


23 


.61803399 


.09289 


.00778 


.00151 


24 


.61803399 


.09289 


.00778 


.00150 


25 


.61803399 


.09289 


.00778 


.00151 



Table 1: Discrepancies for generalized Richtmeyer point sets with iV = 
10 2 , 10 3 , 10 4 points, with for a successive Fibonacci approximations to the golden 
ratio. 
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Figure 2: The function Hq(s) as a function of s. 

of typical integrands. The result of Eq.[22| is useful for large s values, but for 
small values another approach is more efficient. In it, we write z = —u 2 /2, with 
u = a + iv, where a is fixed and v runs over the real axis. This gives us 

G (z) 
H (s) 

H (y; s) 

We can now choose a = yn/s, independently for each m, and so end up with 

i i \/27T / n n \ 7T 2 (2m-l) 2 

#o(s) = £^f (i:\2m-lf -s)e-^^ , (25) 

m>0 S 

which converges very rapidly for small s. This is of course the essence of the 
Poisson summation formula. In Fig. ||] we plot Hq(s) as a function of s. The 
distribution is seen to be considerably skewed, with the maximum around s ~ 1.8 
occuring a good deal before (s) ~ 3.28. 



£F (2m-l;s) , 
) 

dv (a 2 -v 2 + 2iav)e^ a+iv)2 - 7ryia+iv) . (24) 



m>0 
oo 
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The computation of R,2(s;x) runs along the same lines. We may use the 
properties of Bernoulli polynomials to write 



4>{z; x) = 2 \ cos{2imx) = 1 - G (z) cos UVzi^z) , (26) 

n>0 n ~ 2z 7 

with £ = 1 — 2x. Therefore, we may write 
R 2 (s;x) = H {s)-T(s;C) , 

ioo 2 

T(s;0 = — dze~ zs cos (^v^) - 2 ; (27) 

2m J ioo ' (sin(x/2^i)^ 



by closing the contour around the (double) poles, we obtain after some delicate 
but straightforward algebra: 



m>0 



R 2 (s;0 = ]T (-) m ^ 1 m 2 e- sm 12 [l + (m 2 s - 3) cos(27rma;) 

+TTm(2x — 1) sin(27rmx)] . (2<- 



It is easily checked that R 2 vanishes upon integration over either x or s, as it 
should. The form for small s can be obtained along the same lines as above, 
where we may use 



0(-m 2 /2; x)G (z) = G (z) - 2ttV (e™ 5 + e"™ 5 ) me- 2num , (29) 
to arrive at 



m>0 



R 2 (s; x) = H (s)-J2 m l T ( 2m + O + T(2m-O} 

m>0 



T(y) = ^7r 2 t/(7rV-3s)e-— . (30) 



Again, also this form can, in principle, be obtained from Eq.^ by means of the 
Poisson summation formula. In Fig. ^| we show the resulting forms for F 2 (s; x) for 
a number of different values for s. We see that, for small s, this quantity indeed 
peaks at {x} ~ 0. For large s, on the other hand, it is actually negative for small 
x, implying that an anomalously large discrepancy occurs when the points are 
positively correlated, that is, clumped together. Finally, for small s, it is easily 
shown that F 2 (s; 0) ~ tc 2 /s. This gives us an approximate lower bound of ir 2 /N 
on s, which should be compared to Eq.[20|. The bound is quite reasonable, given 
that all our results are obtained in the large- iV limit, where we have assumed 
sN 3> 1. Including more terms ought to enable us to improve on this limit. 

Concerning the low-s approximation, we want to remark the following. We 
can, in principle, also obtain a similar result by a straightforward saddle-point 
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Figure 3: The function i^s; x) for various s. 

approximation, in which the integration over z runs parallel to the imaginary axis, 
with a large negative real part. That this approach must work well is already 
evident from Eq.|29|, in which the leading term is supplemented by exponentially 
small corrections. For large s, a similar saddle-point approximation, where the 
real part of z is close to 1/2, is not so good, since the sub-leading corrections 
are of the same order as the leading term: this yields the correct behaviour with 
s, but with an erreoneous normalization. However, we are mainly interested in 
point-sets with low s, so that we may hope to at least approximate H (s) and 
F 2 (s;a;) by similar saddle-point approximations 0. This however will depend 
sensitively on the mult i- dimensional generalization of a n . 
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5 Conclusions 



In this paper, we have presented, for the first time, expressions for the two-point 
correlation function F 2 (s;xi — £2) for point sets, that are restricted to have a 
certain value s of the discrepancy. These results will be employed to study the 
improvement in integration error for Quasi-Monte Carlo 0. A number of points 
remains to be more fully understood. 

On the one hand, there is the similarity between our form of Hq(s) and the 
Kolmogorov-Smirnov distribution. This is surprising, since they are based on 
quite different discrepancies. 

On the other hand, our 'simple' model lends itself to a generalization to more 
dimensions, but lacks realism. It remains to be seen whether there exists a 
discrepancy-definition, that combines simplicity of results with an easy extension 
to higher dimensions. 
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